! Copyright (c) 2022-2026 Jason Christopherson ! SPDX-License-Identifier: MIT ! ! Permission is hereby granted, free of charge, to any person obtaining a copy ! of this software and associated documentation files (the "Software"), to deal ! in the Software without restriction, including without limitation the rights ! to use, copy, modify, merge, publish, distribute, sublicense, and/or sell ! copies of the Software, and to permit persons to whom the Software is ! furnished to do so, subject to the following conditions: ! ! The Software is provided "as is", without warranty of any kind, express or ! implied, including but not limited to the warranties of merchantability, ! fitness for a particular purpose and noninfringement. module dynamics_vibrations use iso_fortran_env use peaks use ieee_arithmetic implicit none private public :: q_factor public :: estimate_bandwidth public :: logarithmic_decrement public :: damping_from_log_decrement public :: find_free_response_properties public :: rise_time public :: find_settling_amplitude public :: damping_from_fractional_overshoot public :: evaluate_step_response contains ! ------------------------------------------------------------------------------ pure elemental function q_factor(zeta) result(rst) !! Estimates the Q-factor for a vibratory system. The Q-factor is computed !! \(Q = \frac{1}{2 \zeta}\). !! For a lightly damped mode, \(Q\) is also approximately the ratio of the !! resonant frequency to its half-power bandwidth. real(real64), intent(in) :: zeta !! The damping ratio. real(real64) :: rst !! The Q-factor. ! Process rst = 1.0d0 / (2.0d0 * zeta) end function ! ------------------------------------------------------------------------------ pure elemental function estimate_bandwidth(fn, zeta) result(rst) !! Estimates the bandwidth of the resonant mode of a vibratory system. !! The bandwidth is the width of the range of frequencies for which the !! energy is at least half its peak value and is computed as !! \(\Delta f = \frac{f_n}{Q}\). !! Combining this relation with \(Q = 1/(2\zeta)\) gives !! $$ \Delta f = 2\zeta f_n. $$ real(real64), intent(in) :: fn !! The resonant frequency. The units are not important; however, !! the units of the output will be the same as the units of this !! parameter. real(real64), intent(in) :: zeta !! The damping ratio. real(real64) :: rst !! The bandwidth. ! Process rst = fn / q_factor(zeta) end function ! ------------------------------------------------------------------------------ pure elemental function logarithmic_decrement(x1, x2, n) result(rst) !! Computes the logarithmic decrement given the value of two successive !! peaks in the time history of the free vibratory response of the system. !! The logarithmic decrement is calculated as follows. !! !! $$ \delta = \frac{1}{N} \ln \left( \frac{x(t)}{x(t + N T)} \right) = !! \frac{1}{N} \ln \left( \frac{x_1}{x_2} \right) $$ !! For an underdamped SDOF response, !! $$ \delta = \frac{2\pi\zeta}{\sqrt{1-\zeta^2}}. $$ real(real64), intent(in) :: x1 !! The amplitude of the first peak. real(real64), intent(in) :: x2 !! The amplitude of the second peak that occurs N periods after the !! first. integer(int32), intent(in) :: n !! The number of periods of oscillation seperating the two peaks. real(real64) :: rst !! The logarithmic decrement \(\delta\). ! Process rst = (1.0d0 / n) * log(x1 / x2) end function ! ------------------------------------------------------------------------------ pure elemental function damping_from_log_decrement(delta) result(rst) !! Computes the damping ratio from the logarithmic decrement \(\delta\). !! The damping ratio is related to the logarithmic decrement by the !! following relationship. !! !! $$ \zeta = \frac{\delta}{\sqrt{4 \pi^2 + \delta^2}} $$ real(real64), intent(in) :: delta !! The logarithmic decrement. real(real64) :: rst !! The damping ratio. ! Process real(real64), parameter :: pi = 2.0d0 * acos(0.0d0) rst = delta / sqrt(4.0d0 * pi**2 + delta**2) end function ! ------------------------------------------------------------------------------ subroutine find_free_response_properties(t, x, delta, fn, x1, x2, t1, t2, s, n) !! Given a free-response time history, this routine attempts to find the !! logarithmic decrement and resonant frequency of a vibratory system. The !! logarithmic decrement is estimated by finding successive peaks by !! means of peak detection. !! If peaks are separated by \(N\) cycles, the damped frequency estimate is !! $$ f_d = \frac{N}{t_2-t_1}. $$ real(real64), intent(in), dimension(:) :: t !! An N-element array containing the values in time real(real64), intent(in), dimension(:) :: x !! An N-element array containing the response sampled at the time points !! given in t. real(real64), intent(out) :: delta !! The logarithmic decrement estimate. If sufficient peaks cannot be !! located, the routine returns NaN. real(real64), intent(out) :: fn !! The damped resonant frequency in units of Hz, assuming that the !! time values are in seconds. If the time units are not in seconds, !! the units will be cycle/unit time with unit time being the units !! in which t is supplied. If sufficient peaks cannot be located, the !! routine returns NaN. real(real64), intent(out), optional :: x1 !! An optional parameter that, if provided, allows for the routine to !! return the amplitude of the first peak. If sufficient peaks cannot !! be located, the routine returns NaN. real(real64), intent(out), optional :: x2 !! An optional parameter that, if provided, allows for the routine to !! return the amplitude of the second peak. If sufficient peaks cannot !! be located, the routine returns NaN. real(real64), intent(out), optional :: t1 !! An optional parameter that, if provided, allows for the routine to !! return the time at which the first peak was located. If sufficient !! peaks cannot be located, the routine returns NaN. real(real64), intent(out), optional :: t2 !! An optional parameter that, if provided, allows for the routine to !! return the time at which the second peak was located. If sufficient !! peaks cannot be located, the routine returns NaN. real(real64), intent(in), optional :: s !! An optional input that, if provided, allows for control of the !! sensitivity of the peak detection algorithm. The default is 0.1% !! of the peak-peak amplitude of the signal. integer(int32), intent(in), optional :: n !! An optional input that, if provided, determines the number of !! periods to allow between peak selection for the logarithmic !! decrement calculation. The default is 1. ! Local Variables integer(int32) :: np, i1, i2, j2 real(real64) :: xmax, xmin, dx, x1p, x2p, t1p, t2p, nan integer(int32), allocatable, dimension(:) :: maxind, minind real(real64), allocatable, dimension(:) :: maxvals, minvals ! Determine a suitable sensitivity to peak detection if (present(s)) then dx = s else xmax = maxval(x) xmin = minval(x) dx = 1.0d-3 * (xmax - xmin) end if ! Peak Count if (present(n)) then np = n else np = 1 end if ! Additional initialization nan = ieee_value(nan, IEEE_QUIET_NAN) delta = nan fn = nan t1p = nan t2p = nan x1p = nan x2p = nan ! Locate peaks call peak_detect(x, dx, maxind, maxvals, minind, minvals) if (size(maxind) < 2) then ! Return NaN's as we couldn't find enough peaks go to 10 end if i1 = maxind(1) np = min(np, size(maxind) - 1) j2 = np + 1 i2 = maxind(j2) t1p = t(i1) t2p = t(i2) x1p = x(i1) x2p = x(i2) delta = logarithmic_decrement(x1p, x2p, np) fn = np / (t2p - t1p) ! End 10 continue if (present(x1)) x1 = x1p if (present(x2)) x2 = x2p if (present(t1)) t1 = t1p if (present(t2)) t2 = t2p end subroutine ! ------------------------------------------------------------------------------ pure elemental function rise_time(wn, zeta) result(rst) !! Computes the rise time for an underdamped, second-order system. The !! rise time is the time it takes for the system response to go from 0% !! to 100% of its final value and is given by the following relationship. !! !! $$ t_r = \frac{1}{\omega_d} \left( \pi - !! \arctan \frac{\sqrt{1 - zeta^2}}{\zeta} \right) $$ !! where \(\omega_d = \omega_n\sqrt{1-\zeta^2}\) and the formula assumes !! \(\omega_n>0\) and \(0<\zeta<1\). real(real64), intent(in) :: wn !! The resonant frequency of the system, in rad/s. real(real64), intent(in) :: zeta !! The damping ratio of the system. This value must be less than 1 !! as this relationship is only valid for an underdamped system. real(real64) :: rst !! The rise time, in units of seconds. ! Local Variables real(real64) :: arg ! Parameters real(real64), parameter :: pi = 2.0d0 * acos(0.0d0) ! Process arg = sqrt(1.0d0 - zeta**2) rst = (1.0d0 / (wn * arg)) * (pi - atan(arg / zeta)) end function ! ------------------------------------------------------------------------------ pure function find_settling_amplitude(x) result(rst) use fftpack, only : rfft !! Estimates the settling amplitude for a step response. !! The final-value estimate is the zero-frequency Fourier coefficient, !! $$ x_f \approx \frac{1}{N}\sum_{k=0}^{N-1}x_k. $$ real(real64), intent(in), dimension(:) :: x !! The step response of the system. real(real64) :: rst !! The settling amplitude of the step response. ! Local Variables real(real64), allocatable, dimension(:) :: xfft ! Compute the FFT of X and normalize xfft = rfft(x) / size(x) ! We only need the DC component rst = xfft(1) end function ! ------------------------------------------------------------------------------ pure function damping_from_fractional_overshoot(x) result(rst) !! Employs the method of fractional overshoot to estimate the damping ratio !! from the response of a system to a step input. This method is useful !! for cases where the damping ratio is between approximately 0.5 to 0.8. !! In such range, the logarithmic decrement approach becomes less precise. !! !! The fractional overshoot method locates the amplitude of the first !! peak of oscillation (\(x_p\)) and the settling amplitude (\(x_f\)), and !! the estimates the damping ratio as follows. !! !! $$ s = \frac{x_p - x_f}{x_f} $$ !! !! $$ \zeta = \frac{1}{\sqrt{1 + \left( \frac{\pi}{\ln{s}} \right)^2}} $$ !! This follows from \(s=e^{-\zeta\pi/\sqrt{1-\zeta^2}}\) for an !! underdamped step response. real(real64), intent(in), dimension(:) :: x !! The step response of the system. real(real64) :: rst !! The estimated damping ratio. ! Parameters real(real64), parameter :: pi = 2.0d0 * acos(0.0d0) ! Local Variables real(real64) :: xp, xf, s ! Locate the amplitude terms xp = maxval(abs(x)) xf = abs(find_settling_amplitude(x)) s = (xp - xf) / xf ! Compute the damping ratio rst = 1.0d0 / sqrt(1.0d0 + (pi / log(s))**2) end function ! ------------------------------------------------------------------------------ pure elemental function evaluate_step_response(wn, zeta, xs, t) result(rst) !! Evaluates the response of an underdamped single-degree-of-freedom, !! linear system to a step function of amplitude \(X_s\). !! !! The step function response of an underdamped linear SDOF system is given !! as follows. !! !! $$ \ddot{x} + 2 \zeta \omega_n \dot{x} + \omega_n^2 x = \frac{F(t)}{m} $$ !! !! $$ \frac{x(t)}{X_s} = 1 - e^{-\zeta \omega_n t} \left( !! \frac{\zeta \omega_n}{\omega_d} \sin{\omega_d t} + \cos{\omega_d t} !! \right) $$ !! !! where, !! !! $$ \omega_d = \omega_n \sqrt{1 - \zeta^2} $$ !! !! and !! !! $$ X_s = \frac{F}{k} $$ !! Thus the response is the static displacement multiplied by a transient !! term whose envelope decays as \(e^{-\zeta\omega_n t}\). real(real64), intent(in) :: wn !! The resonant frequency, in rad/s. real(real64), intent(in) :: zeta !! The damping ratio. real(real64), intent(in) :: xs !! The amplitude of the step input. real(real64), intent(in) :: t !! The point in time at which to evaluate the response (units = s). real(real64) :: rst !! The step response. ! Local Variables real(real64) :: wd, A ! Process wd = wn * sqrt(1.0d0 - zeta**2) A = zeta * wn / wd rst = xs * (1.0d0 - exp(-zeta * wn * t) * (A * sin(wd * t) + cos(wd * t))) end function ! ------------------------------------------------------------------------------ end module